Exemple simple pour retrouver la variances et la moyenne d’un échantillion gaussien
On utilise ici la version grand dimension et la version grand dimension parallélisé de l’algorithme Metropolis Hastings.
Cet exemple sert à montrer l’utilisation des fonctions MH.
Pour la version parallélisé il est demandé de fournir un écart type sd (possiblement un vecteur d’écart type) mais plus une fonction de transition.
n <- 10000# nombre d'observation
x <- c(3,6) #moyenne et variance
obs <- rnorm(n, x[1], sqrt(x[2])) #tirage dans une loi normale
transi = function(x) rnorm(length(x), x, 0.1) #fonction de transition entre x et xnew
loglik <- function(x, obs) sum( -1/2*log(2*pi*x[2]) -1/(2*x[2])* (obs-x[1])^2 ) #vraisemblence
#Metropolis hasting
x0 <- matrix(c(1,1), nrow = 1)
res <- MH_High_Dim(500, x0, transi, loglik, obs, verbatim = T)
plot(res, nrow = 2)
res_para <- MH_High_Dim_para(500, x0, sd = 0.1 , loglik, obs, verbatim = T, cores = 1)
plot(res_para, nrow = 2)
res %>% as.numeric()
## [1] 3.021748 6.034551
res_para %>% as.numeric()
## [1] 3.013567 5.900558
On peut relancer l’algorithme tout en gardant les informations déjà compiler
res_para2 <- MH_High_Dim_para(100, res_para, sd = 0.1 , loglik, obs, verbatim = T, cores = 1)
plot(res_para2, nrow = 2)
# Exemple d'optimalité entre les 3 fct MH, MH_high_dim et MH_high_dim_para sur le modèle non lineaire mixte
m <- function(t, eta, phi) (phi[1] + eta)/(1+exp((phi[2]-t)/phi[3]))
#=======================================#
param <- list(rho2 = 0.1,
sigma2 = 0.05,
mu = c(5,4,0.5),
omega2 = c(0.5,0.1,0.01))
# mu = 5,
# omega2 = 0.5)
F. <- length(param$omega2) #dimension de phi
#=======================================#
t <- seq(2,6, length.out = 10) #valeur des temps
G <- 40 ; ng = 4
dt <- NLME_obs(G, ng, t, param, m)
Y <- dt$obs
dt %>% ggplot(aes(t, obs, col = gen, group = id)) +
geom_point() + geom_line() +
theme(legend.position = 'null')
#=============================================================================#
Le but ici est de montre la différence de temps de calcul, a noté que certe MH est rapide ne fait pas la même chose que les deux autres fct.
De plus on montre
Phi <- c(- 1/(2*param$sigma2), -1/(2*param$omega2), param$mu/param$omega2)
S <- function(phi)
c(sum((Y - get_obs(m, dt, phi = phi) )^2 ),
apply(phi^2, 2, sum),
apply(phi , 2, sum) )
#=============================================================================#
sd <- 0.02
transi <- function(x) rnorm(length(x), x, sd) #Fonction de transition pour MH_high_dim
transi.noopti <- function(x) rnorm(length(x), x, sd) %>% matrix(nrow = G) #Fonction de transition pour MH
loglik <- function(phi)
{
stat <- S(phi)
sum(Phi*stat)
}
phi0 <- matrix(rnorm(F.*G, c(6,3,1), 0.1), ncol = G) %>% t
n.iter <- 1000
tic()
phi.MH <- MH(n.iter, phi0, transi.noopti, loglik, verbatim = T)
toc()
## 26.79 sec elapsed
tic()
phi.MH2 <- MH_High_Dim(n.iter, phi0 , transi, loglik, verbatim = T)
toc()
## 1063.157 sec elapsed
tic()
phi.MH3 <- MH_High_Dim_para(n.iter, phi0, sd, loglik, verbatim = T, cores = 1)
toc()
## 1057.408 sec elapsed
tic()
phi.MH4 <- MH_High_Dim_para(n.iter, phi0, sd, loglik, verbatim = T, cores = parallel::detectCores() -1)
toc()
## 46.532 sec elapsed
tic()
phi.MH5 <- MH_High_Dim_para_future(n.iter, phi0, sd, loglik, verbatim = T, cores = 1)
toc()
## 1295.729 sec elapsed
tic()
phi.MH6 <- MH_High_Dim_para_future(n.iter, phi0, sd, loglik, verbatim = T, cores = future::availableCores() -1)
toc()
## 91.745 sec elapsed
#=============================================================================#
#Affichage resultat
v <- attr(phi.MH, 'value') %>% lapply(function(x) as.numeric(x)) %>% as.data.frame %>% t %>% as.data.frame()
v %>% mutate(i = 1:nrow(.)) %>% melt(id = 'i') %>% mutate(phi = rep(factor(c(1:F.)), each = nrow(v)*G)) %>%
ggplot(aes(i, value, col = phi, group = variable)) +
geom_line()
plot(phi.MH2, nrow = 2)
plot(phi.MH3, nrow = 2)
plot(phi.MH4, nrow = 2)
plot(phi.MH5, nrow = 2)
plot(phi.MH6, nrow = 2)
#=============================================================================#
n <- 800
m2 <- MH_burnin(phi.MH2, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
summarise(mean = mean(value), sd = sd(value))
m3 <- MH_burnin(phi.MH3, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
summarise(mean = mean(value), sd = sd(value))
m4 <- MH_burnin(phi.MH4, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
summarise(mean = mean(value), sd = sd(value))
m5 <- MH_burnin(phi.MH5, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
summarise(mean = mean(value), sd = sd(value))
m6 <- MH_burnin(phi.MH6, n) %>% attr('value') %>% melt(id = c('id','iter')) %>% group_by(interaction(variable,id)) %>%
summarise(mean = mean(value), sd = sd(value))
mean(abs(m2$mean - m3$mean))
## [1] 0.5902355
mean(abs(m3$mean - m4$mean))
## [1] 0.07936126
mean(abs(m4$mean - m5$mean))
## [1] 0.07876204
mean(abs(m5$mean - m6$mean))
## [1] 0.0958769
mean(abs(m6$mean - m2$mean))
## [1] 0.5659555
mean(abs(m2$sd - m3$sd))
## [1] 0.02676858
mean(abs(m3$sd - m4$sd))
## [1] 0.01533202
mean(abs(m4$sd - m5$sd))
## [1] 0.01184366
mean(abs(m5$sd - m6$sd))
## [1] 0.01640094
mean(abs(m6$sd - m2$sd))
## [1] 0.0326134